Pattern formation in crystal growth under parabolic shear flow II 
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Abstract 

Wavy pattern of ice with a specific wavelength occurs during ice growth from a thin layer of 
undercooled water flowing down the surface of icicles or inclined plane. In the preceding paper [K. 



Ueno, Phys. Rev. E 68, 021603 (2003)], we have found that restoring forces due to gravity and 
surface tension is a factor for stabilization of morphological instability of the solid-liquid interface. 



However, the mechanism for the morphological instability and stability of the solid-liquid interface 
has not been well understood. In the present paper, it is shown that a phase difference between 
fluctuation of the solid-liquid interface and distribution of heat flux at the deformed solid-liquid 
interface, which depends on the magnitude of the restoring forces, is a cause of the instability and 

o 

stability of the interface. This mechanism is completely different from the usual Mullins-Sekerka 
instability due to diffusion and stabilization due to the Gibbs-Thomson effect. 
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FIG. 1: Schematic diagram of vertical cross section of icicle covered with thin layer of flowing water 
and trapped many tiny air bubbles. Mean spacing between ribs is about 1 cm. 

I. INTRODUCTION 

n 

Ripple formation in sand induced by water shear flow [l| and ribs and hollows formation 
on the surface of icicles covered with thin layer of flowing water (see Fig. ^or Fig. 9A in Ref. 
P]) are well known phenomena in nature. The similar wavy pattern as ribs and hollows on 
icicles in nature can be experimentally produced during ice growth by continuously supplying 
a proper water Q ml/h on an inclined plane with width I and at angle 8, set in cold room 
below 0° C sketched in Fig. and it is found that the mean wavelength of the wavy pattern 
of ice is given by 0.83/(sin #)°- 6 ~ - 9 cm In Fig. |2l the shaded regions with uniform spacing 
are protruded part of the wavy pattern of ice. Indeed, the spacing of the wavy pattern at 



6 = Ti / 2 is nearly equal to t 



In the previous works 
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le mean spacing between ribs on the surface of icicles in nature. 



a morphological instability of the solid-liquid interface during 
a crystal growth with mean velocity V from an undercooled thin liquid flowing down an 
inclined plane under the action of gravity as shown in Fig. El was investigated. We restricted 
ourselves to two dimensions in a vertical plane (x, y) , and for simplicity we assumed that 
the region of the crystal is semi-infinite. The x axis is parallel to the inclined plane and the 



y axis is normal to it. The parabolic shear flow 



is parallel to the x axis and is bounded on one side by a liquid-air surface which is exposed 
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FIG. 2: Schematic diagram of ice growth from thin layer of undercooled water flowing down under 
the action of gravity on inclined plane with width I cm and at angle 8, which is set in a cold room 
below 0° C. Water supply rate is Q ml/h. Shaded regions with uniform spacing are nearly periodic 
wavy pattern of ice. 

by cold air below 0° C. Here uq is the velocity at the free surface and ho is the mean thickness 
of the liquid, which can be expressed as 



where g is the gravitational acceleration and v is the kinematic viscosity. 

We discussed some differences between the dispersion relation for the fluctuation of the 
solid-liquid interface in our model and that in the Ogawa and Furukawa's model ( hereafter, 
we refer to thier model as O-F model ). Our amplification rate a r and phase velocity v p are 
given by jj] 




(2) 



V -|a(//Pe)+/x{36-§a(/xPe)} 
h 36 + a 2 



+ njj, 



j^a(/zPe) — a 2 + fi {36 — ^a(/zPe)} 
36 + a 2 



(3) 



V -\a 2 {fiPe) + /i{6a + 9(/iPe)} 
ji 36 + a 2 



+ n\i 



6a - ^a 2 (/xPe) + fj, {6a + f QuPe) } 
36 + a 2 



(4) 
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FIG. 3: Schematic diagram of vertical plane (x, y) of inclined plane at angle 0. 



where \i = kho, Pe = uqIio/ki is the Peclet number, n = K s /Ki, k, ki, K\ and K s being the 
wave number, the thermal diffusivity of the liquid, the thermal conductivity of the liquid 
and solid, respectively. 

a = 2 cot 9h k + a 2 h Q k 3 (5) 

is restoring forces due to gravity and surface tension acting on the liquid-air surface, a = 
sj^ff[fiig sin 9) being the capillary constant associated with the surface tension 7 of the 
liquid-air surface. While, by using our notations, the O-F model gives \^ 
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(//Pe) 



a r = Vk 10080 V ~ 2: (6) 

{l-lilo(^} 2 +{^Pe} 2 ' 

v p = V 3. 7 

{l-lilo(^} 2 + {^Pe} 2 

A critical differece of these dispersion relations is that our result includes the restoring 
forces a, while the result of O-F model does not include it. The difference of a r results 
in different dependence of wavelength A max , at which <x r takes a maximum value, on 9. 
Our result is better agreement with A mea n obtained by experiment f.3] than the result of 
O-F model (see Fig. 4 in Ref. ^\)- The difference of v p results in different direction of 
migration of the solid-liquid interface. Our results predict that it moves upward with about 



0.6V. Indeed, our prediction is consistent with the observation that many tiny air bubbles 
trapped in just upstream region of any protruded part migrate in the upward direction 
during growth as shown in Fig. ^ On the other hand, the result of O-F model predicts 
that it moves downward with about 0.5V. Since our calculations based on a linear stability 
analysis and the O-F's calculations are different, the cause leading to these different results 
was not clarified in the previous paper. In addition to the presence of a or not, it was only 
suggested that these differences may be due to the difference of boundary conditions of the 
temperature at the solid-liquid interface and that of the liquid-air surface. In Sec. II, we 
confirm it by deriving O-F model from our formulation. 

According to the result of O-F model, the instability of the solid-liquid interface occures 
by the Laplace instability due to the thermal diffusion into the air, and its instability is 
suppressed by the effect of fluid flow, which makes the temperature distribution in the thin 
water layer uniform [5]. However, this qualitative interpretation does not enable to explain 
the migration of the interface, and the role of fluid flow on the temperature field in the 
liquid is not quantitatively clear. In our previous paper, although it was suggested that the 
restoring forces are indeed an important factor for stabilization of morphological instability 
of the solid-liquid interface (see Fig. 2 in Ref. the mechanism was not well understood. 
In Sec. Ill, therefore, we clarify the morphological instability and stability mechanism of the 
solid-liquid interface, and we present a physical mechanism for migration of the interface in 
the upstream direction. Conclusion is given in Sec. IV. 

II. DIFFERENCE IN BOUNDARY CONDITIONS 

The general solution for the amplitude of perturbed temperature in the liquid is 



9l (z) = BifaW + B 2 cf) 2 {z) + ifiPe / {M*)<hW) - <h{z)<hW)} fiz'WGta, (8) 



where 4>i(z) and <f>2(z) are given by Eqs. (81) and (82) in Ref. |4(, Gi is the unperturbed 
temperature gradient in the liquid, and z — 1 — y/ho. 

The first difference of our model and O-F model is the form of the amplitude / of per- 
turbed stream function in Eq. (JHJ). The solution of / was determined under the same 
hydrodynamic boundary conditions in the two models, but the difference of the order esti- 
mate of the restoring force a led to different solutions. Our result in the long wavelength 
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approximation /i<l gives 

f(z) = (— 6 + iaz + 6z 2 — iaz 3 ) 

6 — la 

fiRea , _ 2 

210(6 -ia) 21 V 7 

+ (210 - Ilia)/ + (-42 + 7ia)z 5 + (6 - iar)* 7 } , (9) 

where Re = uoho/u is the Reynolds number. While the result of O-F model's gives 

f(z) =z 2 -l. (10) 

From the kinematic condition at the liquid-air surface, the reletion between the amplitude 
(k of the perturbed solid-liquid interface and the amplitude of the perturbed liquid-air 
surface was given by = — f\ z= o(k 4]. Equation © includes a which makes the amplitude 
and phase of the liquid-air surface change, it deforms depending on the wavelength of the 
fluctuation of solid-liquid interface. While, Eq. (jlOj) does not include it. Then, the liquid-air 
surface fluctuates with the same amplitude as the solid-liquid interface and phase shift of 
each interface does not occur. 

The second difference between our model and O-F model is originated from the difference 
of B\ and B2 in Eq. (|Sji . which depend on the choice of the boundary conditions of the 
temperature at the solid-liquid interface or the liquid-air surface. The two differences lead to 
a critical difference of the dispersion relations for the fluctuation of the solid-liquid interface 
and of the corresponding mechanism of instability and stability of it. The difference of the 
mechanism will be discussed in detail in the next section. Here we describe the essential 
difference of the thermodynamic boundary conditions in the two models, and we correct the 
solution Eq. (86) of unperturbed temperature of the air obtained in Ref. 

The first essential difference is the continuity of the temperature at the perturbed solid- 
liquid interface C(tj x ) = Cfc^Pt^ + ikx], where a = a r + icr, and t is time. In our model, 



Ti\ y =£ = T s \ y=i = T m + AT, (11) 

where T m is the equilibrium melting temperature. Here we assume that a deviation from the 
equilibrium melting temperature, AT, is of order (k and this corresponds to G(k)( in the 
previous paper 4] . If we regard AT as the melting temperature depression due to the Gibbs- 
Thomson effect ; 8] , we can neglect it as far as we are concerned with the wavelength of the 



wavy pattern observed on the surface of icicles or the inclined plane. Then, the temperature 
at the solid-liquid interface in a pure substance must be the equilibrium melting temperature. 
Therefore, in the O-F model, 

Ti\ y= c_ = T 8 \y=£ = T m . (12) 

However, we suggested that there can be a deviation from the equilibrium melting tempera- 
ture, which can not be determined a priori and is determined after we determine the solution 
for the perturbed temperature in the liquid 0. 

The second essential difference is the continuity of the temperature at the liquid-air 
surface x) = h + expfcrt + ikx]. In our model, 

Ti\ y= £ = T a \ y= £ = Tic, (13) 

where 7] a is a temperature at the liquid-air surface. While, in the O-Fmodel, 

T t \ y=( = T a \ yH . (14) 

Equation (|T3*j) indicates that the temperature of the liquid-air surface remains a constant 
Ti a after deformation of the liquid-air surface. While, Eq. ()14|) shows that the temperature 
at the deformed liquid-air surface is not necessary to remain a constant value. 

From the heat conservations Eq. (94) at the unperturbed liquid-air surface and Eq. (101) 
at the unperturbed solid-liquid interface in Ref. we obtain TJ = T 00 + Lj 'C pa , T^, L and 
C pa being the ambient air temperature, the latent heat per unit volume and the specific heat 
of the air at constant pressure, respectively. In this case, nothing determines the growth 
velocity V, which may take any value. This problem also occur in the absence of flow 0. 
There is a more serious problem. The value of LI Cr,„ is about 2.54 x 10 3 K. Then, for the real 

nn 

values of observed [2, |3J|, the value of 7] a becomes larger than the melting temperature. 
This problem is originated from the unperturbed solution (86) of the temperature T a of the 
air in Ref. [^J. Here we briefly improve the unperturbed solution. 

Under the boundary conditions T a = T ia at y — h and T a = at y = h + l a , we 
assume that by neglecting the term V in Eq. (85) in Ref. the approximate solution of 
the unperturbed part near the liquid-air surface is given by 

T a (y) = T la -G a (y-h ), (15) 

where G a = (Ti a — T^/la is the unperturbed temperature gradient in the air at y = h , l a 
being a length of thermal diffusion layer ahead of the liquid-air surface. We note that this 



l a is not l a = n a /V in the previous paper |4|. From the heat conservation K[G\ = K a G a at 
the unperturbed liquid-air surface, where K a is the thermal conductivity of the air, Ti a is 
obtained as 

rp | Kg hp rp 

™ ~ If, I ™ , 

T - = • ( 16 ) 
^ K t l a 

If the values of T^, /i and l a are given, the value of T la in Eq. (jT3Jl is determined. Sub- 
stituting Eq. into the heat conservation LV = K\Gi at the unperturbed solid-liquid 
interface yields 

Since Ki/K a 3> 1, Eq. (fTTj) can be approximated as 

K T —T 

,/ss -^^' (18) 

We note that V" does not depend on Hq which varies with Q. Therefore, V is not affected 
by change of Q. This agrees with observations that as increasing water supply rate growth 
velocity of diameter of icicles is almost constant |2|. In the previous paper, we determined 
the value of T\ a from only Eq. (115) in Ref. |4j using the actual observed value of V. 
However, the expression (115) is not appropriate for determinig V because it changes with 
h - 

Under the boundary conditions (jllj) and ([13)1 . we obtain j4j 

Bi = -fl^QiCk, B 2 = fiB u (19) 

where f\ z= o is given by Eq. (125) in Ref. 0. Then, the dispersion relation for the fluctuation 
of the solid-liquid interface becomes 

V (dH t 



a 



ho { dz 
where 



+ nfi(H l \ z=l -l)\, (20) 

z=l 



Hi{z) = -f\ g=0 {<l> 1 (z)+^h(z)} + ifJiPe [ {faWM*) ~ Mz)<hW) WW- (21) 

The real and imaginary part of Eq. (|2T)|) give Eqs. (JHJ) and (jlj) by approximating Eqs. (97), 
(119) and (120), respectively 0. In the previous paper, nevertheless we had inappropriate 
base state of the temperature of the air, we obtained reasonable results compatible with 
experiments and observations. This is because the change of V in Eq. (|TH|) by affects 
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FIG. 4: The amplification rate a r vs the wave number k for V = 10 6 m/s, Q = 160 ml/h, and 
9 = it/2. Solid line: Re[Eq. (JUJ)]. Dashed line: Re[Eq. ©]. Dotted line: Re[Eq. (a = 0). 

the magnitude of <r r but does not make change the characteristic wavelength of the wavy 
pattern determind from the maximum point of a r . 

On the other hand, under the boundary conditions (fT2|) and ()14jl . we obtain a different 



Bl = ~~ — ^u+W*iu "' e,a ' (22) 



B<1 = ri — , — ktit\ ( 23 ) 

01 1 2=1 + ^4>2\z=l 

where 

= f {0 2 (*)0i CO - M*)M»?)} fVW- (24) 

Then, the dispersion relation becomes 

V dHi 



a 



ho dz 
where 



z=l 



(25) 



Hi(z) = — ^ — [{1 - ifiPeI\ z=1 } 4>i(z) + nf\ z=0 {fa\ z =i<t>i(z) ~ ^>i\»=ih( z )}] + ifiPeI(z). 

0l|z=l 

(26) 

In Eq. (126)1 . we have omitted the term nK a /Ki <C 1. In particular, when puttinig a = in 
/ in Eq. (J2EJ), we recover Eqs. © and (J7J) from the real and imaginary part of Eq. ([25)1 . 
respectively. 
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FIG. 5: Phase velocity v p = —(Ji/k vs the wave number k for V = 10 6 m/s, Q = 160 ml/h, and 
6 = it/2. Solid line: -Im[Eq. pUJl ]/k. Dashed line: -Im[Eq. l f 2H )l ]/k. Dotted line: -Im[Eq. l|351l ]/k 
(a = 0). 

In the absence of flow, if we regard AT in Eq. (jllj) as the Gibbs-Thomson effect, and by 
replacing the amplitude relation = — f\ z=0 ( k |4( with 



we obtain 



£ k = exp(-/i) (l - d ^k 2 ^) Cfe, 
Hi(z) = exp[-//(l - *)] (l - rf ^ 2 ) • 



(27) 



(28) 



Then, Eq. fOU|) obtained from the boundary conditions (fTT|) and (|T^|) reduces to the disper- 
sion relation in the Mullins-Sekerka theory ^]: 



Vk 



1 - d ^{l + n)k 2 



(29) 



and v p = 0, where d Q = T m TC p /L 2 is the capillary length, T being the solid-liquid interface 
tension jfj. We note that we can not recover the dispersion relation in the Mullins-Sekerka 
theory from the boundary conditions ()12j) and ()14J) even if we add the Gibbs-Thomson effect 
to Eq. (DP). 

The solid lines, dashed lines and dotted lines (a = 0) in Figs. E| and El represent a r and 
Vp = -Gi/k obtained from Eqs. flUJ) and J2SD with the use of Eqs. (97), (119) and (120), 
respectively. Here Re and Im denote the real and imaginary part of arguments, respectively, 
and note that Re is not the Reynolds number in the following discussion. It is found from 
the dashed line in Fig. |U that a r is always positive, therefore, it is not impossible to get 
the characteristic wavelength observed on the surface of icicles or the inclined plane. The 
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dotted line in Fig. 0]has a maximum point of a r at a specific wave number. However, since 
this result of 0-F model does not consider the effect of restoring forces due to gravity and 
surface tension on the liquid-air surface, deviation from the experiment is large as shown in 



the closed triangles in Fig. 4 in Ref. 



Q. 



As shown in the solid line in Fig. the direction 
of phase velocity in our model is negative in our interest wave number region. On the other 
hand, as shown in the dashed and dotted lines in Fig. the direction of phase velocity in 
the models different from ours is positive in our interest wave number region. There is no 
evidence to support this prediction. 

From these considerations, in order to explain the experiments and observations, and in 
the absence of flow, to recover the dispersion relation in the Mullins-Sekerka theory, the 
boundary conditions and (J13|) seem to be most appropriate. 
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The deviation from the equilibrium melting temperature in our model is given by 

AT=(^| 2 =i-l)^C, (30) 

where Hi is given by Eq. (}2"T]) . In the absence of flow, instead of Eq. (}2"T]) . if we apply Eq. 
()28|) to Eq. ()3CH). we recover the Gibbs-Thomson effect, AT = —dolik 2 Gi(, where li = K\jV ' . 
We note that as k — > 0, the solid-liquid interface is flat and both Hi\ z= i in Eqs. (12 1|) and (|2*H|) 
approach 1, therefore, AT vanishes. This indicates that the deviation from the equilibrium 
melting temperature in our systems is not induced by only the effect of shear flow Eq. 

We note the difference of the characteristic length scale in our problem with flow and in 
the Mullins-Sekerka theory. The capillary length do associated with the solid-liquid interface 
tension is a microscopic length of order angstroms, while li is usually macroscopic. Therefore, 
the Gibbs-Thomson effect acts effectively on the micrometer scale On the other hand, 
the restoring forces in Eq. (0) include the capillary constant a assosiated with the surface 
tension of the liquid-air surface, which is 3.9 mm for water at 9 = n/2, and the thickness ho 
of the liquid, which is about 10~ 4 m. Then, the effect of restoring forces is more effective 
for longer wavelength fluctuation compared to the length scale the Gibbs-Thomson effect is 
effective. We can neglect AT due to the Gibbs-Thomson effect but not neglect it due to the 
restoring forces. Indeed, the effect of AT is reflected in the second term of Eq. (}20|) . 

Local equilibrium thermodynamics assumes that the equations of state retain the same 
form out of equilibrium as in equilibrium, but with a local meaning jjj. Then, the flow does 
not change the equations of state and the equilibrium melting temperature is determined 
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from equality of the chemical potentials of the crystal and liquid. This local equilibrium 
hypothesis is implicitly used in the O-F model and leads to Eq. (jl2J) . However, if we impose 
the boundary condition (fT2|) in our systems, we can not get desired results compatible with 
experiments and observations. The restoring forces a have a purely hydrodynamic origin, 
but there is an interplay of both thermodynamic and hydrodynamic effect in AT determined 
from Eq. (j2~Tj) . which is reflected in the term aPe in Hi. 

The thermodynamics of fluids under shear flow is an active and very challenging topic in 
modern nonequilibrium thermodynamics and statistical mechanics jjj. This is a field with 
many open questions. A decisive step in the thermodynamic understanding of AT is to 
formulate a free energy or a chemical potential depending explicitly on the characteristics 
of the flow. However, it is reported that the coexistence of the crystal and liquid under 
shear flow cannot be accounted for by invoking a nonequilibrium analogue of the chemical 
potential Q|. We have seen that if we impose the equilibrium melting temperature at 
the solid-liquid interface, we can not get consistent results with actual experiments and 
observations. Therefore, the approach by the local equilibrium hypothesis may also be 
insufficient to deal with our systems. Whether the coexistence of the crystal and liquid 
under shear flow in the present case can be accounted for by introducing a nonequilibrium 
analogue of the chemical potential is left for future studies. 

III. MECHANISM OF INSTABILITY AND STABILITY OF THE SOLID-LIQUID 
INTERFACE 

In the preceding section, we have studied some models which differ in the form of bound- 
ary conditions. In this section, we present in detail the mechanism of instability and stability 
of the solid-liquid interface for our model. For other models, we give a brief discussion. 

As shown in Fig. |f3 we consider a small perturbation of the solid-liquid interface with 
Ck/ho = 0.1 at time t: 



Using £fc = — f\ z=0 ( k |4(, the corresponding perturbation of the liquid-air surface is 



Im — = sm[k(x — v p t)]— exp(a r t). 
_h j h 



(31) 



Im — = Im [-f\ z=0 exp[ik(x - v p t)]) — exp(oyt). 



(32) 
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Next we define the perturbation of heat flux into the liquid qi and solid q s at the perturbed 
solid-liquid interface and the perturbation of heat flux into the air q a at the perturbed liquid- 
air surface as follows, respectively, 

dT( 



Qi 



Im 
Im 



-Ki 



dH, 



dz 



dy 



z=l 



y=C. 

exp[ik(x — v p t)} 



h 



(33) 



q s = Im 



dy 



-lm[(Hi\ z= i - 1) exp[ik(x - v p t)]}njjKiGi— exp(a r t), (34) 

ho 



q a = Im 



dT' 




K 

dy 





Im [-f\ z=0 exp[ik(x 



h 



(35) 



where we use Eq. (J2T|) for Hi(z). We note the direction of heat flow. If q\ > or q s < 0, the 
latent heat is released away from the solid-liquid interface into each phase. Conveniently, 
the distribution of q\ — q s represented by the bottom dashed line and q a represented by the 
up dashed line are superimposed on Fig. El with magnification of 0.1 to see phase difference 
between the fluctuation of the solid-liquid interface, liquid-air surface and distribution of 
heat flux at the respective interfaces. 

Figure. El (a) shows the configurations at a wave number in the unstable region a r > of 
the solid line in Fig. |3J The heat flux q a at the liquid-air surface is large at any protruded 
part of surface pointing into the air, at which the temperature gradient increases so that 
heat transfer by thermal diffusion into the air is more effective. Since the value of a is small 
in such low wave number region, that is, the effect of restoring forces on the liquid-air surface 
is small, the liquid-air surface fluctuates with almost the same amplitude as the solid-liquid 
interface, and the phase difference between Im[£/Ao] and Im[£/7io] is negligible. Therefore, 
this seems to result in faster cooling and hence freezing to promote at the protruded part 
of the solid-liquid interface. This picture of destabilization appears to be the same as the 
Mullins-Sekerka instability |8|. However, we note that the maximum point of heat flux 
qi — q s is shifted to the upstream direction by against the solid-liquid interface. This 
indicates that the interface grows faster in just the upstream region of any protruded part, 
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FIG. 6: Schematic illustration of the fluctuation of the solid- liquid interface Im[£//io] (bottom 
thick solid lines), the liquid-air surface 1 + Im[£//io] ( U P thick solid lines), heat flux qi — q s (bottom 
dashed lines) at Im[£//io], and heat flux q a (up dashed lines) at 1 + Im[£//io] for (a) k = 634/m, 
(b) k = 953/m, and (c) k = 1200/m. (f> is phase shift of heat flux qi — q s against the solid-liquid 
interface. 
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in which q\ — q s is large compared to the mean heat flux K\Gi. On the other hand, in the 
downstream region of any protruded part, the interface tends to melt back because q\ — q s 
is small compared to K\Gi. Therefore, the solid-liquid interface not only grows unstably 
but also moves in the upstream direction, which is consistent with the direction of phase 
velocity of the interface predicted by Eq. P|). It also support the observation that many tiny 
air bubbles dissolved in the thin flowing water are trapped in just upstream region of any 
protruded part on a growing icicle and its region migrates in the upward direction during 
growth as shown in the dotted regions in Fig. ^ We can not explain this observation by 
usual Mullins-Sekerka instability or Laplace instability [^| due to diffusion. 

Figure El (b) shows the configurations at the wave number at the neutral stability point 
o r = of the solid line in Fig. 0] Shift of the maximum point of heat flux qi — q s is larger 
than that in Fig. El (a). Figure El (c) shows the configurations at a wave number in the stable 
region <j r < of the solid line in Fig. |3J Since the value of a increases as increasing the wave 
number, the surface Im[£//i ] and heat flux q a is slightly shifted to the upstream direction, 
and the phase difference between Im[£//io] and q\ — q s becomes larger than that in Fig. El 
(b). Any protruded part of the perturbed solid-liquid interface melts back because q\ — q s is 
small, and any depression part of the interface grows because qi — q s is large. Therefore, the 
flatness of the solid-liquid interface is restored, that is, the solid-liquid interface is stabilized 
by large phase shift of distribution of heat flux by large restoring forces. This stabilizing 
mechanism on such long length scales shown in Fig. El is different from the Gibbs-Thomson 
effect. 

If we apply the picture of the Mullins-Sekerka instability to the bottom thick solid lines 
in Fig. El the deformed isotherms get closer to each other ahead of the bump of the solid- 
liquid interface. The temperature gradient, and therefore the heat flux, increases, which 
increases the rate of production of latent heat. Therefore, the bump must get amplified jsj]. 
However, we note that this picuture is true only in the absence of flow. In the Mullins- 
Sekerka instability, it is diffusion which destabilizes the planar front. In the presence of 
flow, the perturbed temperature field in the liquid is affected by the flow field, which varies 
depending on the magnitude of the restoring forces acting on the liquid-air surface. We can 
not determine a priori where the temperature gradient or heat flux is large until we solve 
the equation of the temperature field for a given boundary conditions. 

This mechanism of instability and stability of the solid-liquid interface and its movement 
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to the upstream direction discussed above can be explained more quantitatively as follows. 
The perturbed part of heat conservation equation (18) in Ref. ^] is 

-K 9T! 



dt s dy 



dy 



(36) 



Taking imaginary part of Eq. (}36|) . it becomes 



LIm[crexp(crt + ikx)](k = qi — q s , 



(37) 



where qi and q s are given by Eqs. (J33)) and (J34j) . The imaginary part of the left hand side 
of Eq. (|3*Tj) can be written as 



Im[crexp(cr£ + ikx)] = \cr\ exp(o>t) sin [k(x — v p t) 



where \a\ 



'at 



erf 



and 



Or 



cr cos ( 



(Ti 



a sm i 



(38) 



(39) 



Figure El (a)- (c) show that <fi < 0. Noting that v p = —<7i/k, from the second equation of Eq. 
(jnHJ), when < 0, (jj is positive, therefore, v p < 0. From the first equation of Eq. (J3~9"|) . we 
find the sign of <x r and the corresponding figures as follows: 



>o (-!< 



< o ) Fig. EN 



Or < 



-f ) Fig. EKb) 



(40) 



< ( -7T < 6< -\ 



Fig. Elc) 



Equation indicates that unstable, neutral, and stable regions of the solid line in Fig. H] 
are completely consistent with Fig. E](a)-(c), and that the direction of phase velocity is the 
same as the solid line in Fig. El if we restrict ourselves to the wavelength region observed on 
the surface of icicles or the inclined plane. 

Likewise, if we apply Eq. (ffijjl for Hi(z) to Eqs. and (JHIJ), > 0, therefore, 

we obtain v p > 0. This is consistent with the direction shown in the dashed and dotted 
lines in Fig. El According to the O-F model, the stability of the solid-liquid interface 
is due to uniformalization of the temperature distribution along the layer by fluid flow 0. 
However, the fluid flow never make uniform the temperature distribution. If we give a correct 
interpretation for the stabilization of the solid-liquid interface in the O-F model, which is 
essentially the same as that explained in our model. In the stable region, fluctuation of the 
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solid-liquid interface and distribution of heat flux tend to be out of phase. At the protruded 
part, heat flux is small, while at the depression part, heat flux is large, therefore the flatness 
of the interface is restored. As a result, the dotted line in Fig. HJis obtained. For the dashed 
line in Fig. out of phase between fluctuation of the solid-liquid interface and distribution 
of heat flux never occur. Therefore, this case is always unstable. 

IV. CONCLUSION 

We have provided a physical interpretation for the morphological instability and stability 
of the solid-liquid interface occurring during a crystal growth from an undercooled thin 
parabolic shear flow of water on the surface of icicles or the inclined plane. The wavy 
pattern with a characteristic wavelength are observed on longer length scales compared to 
the one determined by the competition of the Mullins-Sekerka instability due to diffusion 
and stabilization due to the Gibbs-Thomson effect. We have found that phase difference 
between fluctuation of the solid-liquid interface and distribution of heat flux at the deformed 
solid-liquid interface, whose difference depends on the magnitude of restoring forces due to 
gravity and surface tension, is the cause for destabilization or stabilization of the interface, 
and that the direction of phase shift of the distribution of heat flux against the solid-liquid 
interface is related to the direction of migration of the solid-liquid interface. 
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